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ABSTRACT 

A  sea  swell  and  surf  model  is  implemented,  tested  and 
evaluated  on  a  micro  computer  (HP-9845B) .   The  two-dimensional 
model  includes  spectral  wind-wave  generation  in  open  water, 
and  shallow  water  wave  transformation  over  irregular 
topography.   The  model  predicts  surf  zone  width,  breaker 
lines  and  types  of  breakers.   Using  change  in  momentum  flux 
of  gravity  waves  (radiation  stress)  as  forcing,  the  model 
predicts  current  velocities  within  the  surf  zone.   The  model 
is  evaluated  for  the  conditions  over  constant  depths  and 
uniform  sloping  beaches.   The  numerical  results  are  checked 
against  accepted  theory  and  field  observations.   The  model  is 
found  to  overbuild  wind  generated  sea  heights  for  a  30  kts 
wind  but  to  give  expected  wave  heights  for  a  15  kts  wind. 
The  model  results  compare  well  with  observed  nearshore  wave 
heights  but  give  poor  location  of  breaking  waves.   The 
model's  nearshore  current  calculation  is  found  unsuitable  for 
the  HP-9  845B  due  to  computational  instability  and  time 
requirements.   The  sources  of  the  model's  problems  are 
identified,  and  recommendations  are  made  for  future  improve- 
ments of  the  model. 
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I.   INTRODUCTION 

A.   MOTIVATION 

Wave  and  surf  forecasting  is  a  primary  concern  of  the 
U.S.  Navy  in  executing  amphibious  operations.   The  wave 
environment  affects  the  larger  vessels  in  deep  water  and 
landing  craft  approaching  the  beach.   As  waves  propagate  from 
deep  to  shallow  water,  thev  undergo  shoali  ;,  refraction,  and 
eventually  breaking,  while  forming  a  surf  zone.   The  height 
and  type  of  breakers  in  the  surf  zone  determines  the  difficulty 
that  landing  craft  will  have  in  reaching  the  beach.   The 
breaking  waves  can  also  generate  nearshore  currents;  these 
shear  currents  potentially  can  turn  a  craft  and  cause  it  to 
be  broached  and  swamped  by  breaking  waves.   In  addition  to 
amphibious  units,  coastal  structures  and  pierside  ships  are 
subject  to  forces  by  waves  in  shallow  water  regions. 

The  forecasting  of  wave  conditions  can  be  divided  into 
two  regimes:   the  nearshore,  including  the  surf  zone,  and  an 
offshore  area.   There  are  a  number  of  numerical  models  avail- 
able which  predict  waves,  shear  currents  and  surf  conditions 
separately.   Propagating  wave  energy  from  deep  water  to  the 
high  water  line  on  the  beach  requires  a  combination  of  comple- 
mentary models.   The  objective  of  this  thesis  is  to  convert  a 
state-of-the-art  computer  program  for  use  on  a  Hewlett-Packard 
9845B  mini  computer,  to  test  the  model  against  field  observations 


and  to  revise  where  necessary.   The  program  will  be  submitted 
as  a  prototype  wave  and  surf  model  for  shipboard  use  in  the 
Fleet. 

B.   BACKGROUND  OF  WAVE  AND  SURF  PREDICTION 

The  earliest  models  of  wave  generation  and  shallow  water 
transformations  were  for  monochromatic  regular  waves.   Refrac- 
tion was  based  on  the  wave  ray  method  (Arthur,  19  49) ,  where 
ray  paths  were  calculated  for  individual  waves  with  slowly 
varying  wave  speeds  along  the  wave  front.   This  method  is 
analogous  to  ray  path  prediction  used  in  optics  for  light 
rays.   With  the  development  of  computers  and  a  better  under- 
standing of  the  physical  processes,  it  became  possible  to 
consider  the  entire  wave  spectrum  in  generation  and  shallow 
water  transformation.   Karlsson  (1969)  introduced  a  direction 
variable  to  study  wave  spectral  transformation  for  straight 
and  parallel  contours.   Collins  (1972)  included  bottom  friction 
and  evaluated  the  validity  of  one  dimensional  models.   Noda, 
et  al.  (1974)  used  a  relaxation  finite-difference  scheme  to 
solve  the  stationary  wave  spectral  transformation  in  shallow 
water  with  the  inclusion  of  bottom  dissipation  and  effects  of 
local  wind  generation.   Shiau  and  Wang  (1977)  studied  wave 
transformation  in  shallow  water  and  compared  the  results  with 
field  observations.   Wang  and  Yang  (1981)  applied  their  model 
with  the  inclusion  of  bottom  friction.   Recently,  Chen  and 
Wang  (1983)  improved  the  model  by  considering  non-stationary 
waves  and  interactions  with  currents. 


C.  CURRENT  STATUS  OF  MODELING  TECHNIQUE 

The  current  state-of-the-art  for  numerical  modeling  can 
be  characterized  by  two  numerical  methods,  finite  difference 
and  finite  elements.   Both  methods  have  been  applied  to  wave 
and  surf  problems.   A  finite  element  model,  such  as  by  Wu  and 
Liu  (1985)  ,  economically  handles  localized  irregular  topography 
and  has  good  properties  for  nonlinear  convergence.   However, 
finite  element  schemes  are  highly  complex  and  require  a  sophis- 
ticated understanding  of  the  model  chara  teristics  to  choose 
the  best  grid  system.   Finite  difference  models  (Noda,  19  74; 
Shiau  and  Wang,  1977)  use  rectangular  grids  in  general,  and 
the  numerical  analysis  has  been  widely  studied.   The  numerical 
procedure  of  finite  difference  models  is  self-explained.   There- 
fore, a  finite  difference  scheme  is  less  demanding  on  the 
user's  numerical  skills  and  from  an  operational  point  of  view, 
is  more  practical. 

D.  MODEL  OVERVIEW 

The  model  described  below  was  developed  by  Wang  and  Chen 
(1983).   It  is  a  finite  difference  model  which  considers  a 
wind  generated  spectrum,  spectral  and  swell  wave  transforma- 
tions over  moderately  irregular  bottom  topography,  and 
dissipation  by  breaking  waves  and  bottom  friction.   The  model 
calculates  a  wind  generated  and  deep  water  wave  spectrum  in 
the  ocean  and  predicts  wave  and  surf  conditions  on  the  beach. 
There  is  also  an  option  of  direct  input  of  a  predetermined 


wave  into  a  surf  area.   In  the  nearshore  computations,  th  -re 
are  options  for  computing  surf  conditions  and/or  nearshore 
circulation.   The  surf  and  nearshore  circulation  output  are 
the  basis  for  predicting  amphibious  landing  conditions.   Model 
input  includes  wind  direction  and  speed,  initial  ocean  swell 
and  dominating  direction,  bottom  type,  and  water  depths.   All 
the  inputs  are  entered  for  the  deep  water  calculation  and 
again  serve  as  appropriate  data  for  the  shallow  water. 

E.   OBJECTIVES 

The  objective  of  this  thesis  is  to  compile  the  numerical 
model  by  Wang  and  Chen  (19  83)  in  a  user  friendly  version  and 
test  it  on  the  HP-9  845B  mini  computer.   The  program  should  be 
self-explanatory  with  internal  error  checks  and  simple  correc- 
tion steps.   The  model  will  be  tested  against  observations  on 
several  beaches  to  establish  accuracy  and  range  of  application 
The  physical  and  numerical  limitations  of  the  model  will  be 
documented,  explained  and  where  possible  improved.   It  should 
be  a  useful  tool  for  the  Shipboard  Aerographer  Mates. 
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II.   THEORY  AND  STRUCTURE  OF  THE  MODEL 

A.  THE  SEA  SWELL  AND  SURE  PROGRAM  (SSSP) 

The  sea,  swell  and  surf  program  (SSSP)  was  developed  by 
Wang  and  Chen  (1983)  and  is  a  combination  of  several  comple- 
mentary models.   The  sea  and  swell  for  open  sea  follows  the 
models  by  Shiau  and  Wang  (197  7)  and  Chen  and  Wang  (19  82) .   The 
nearshore  current  calculation  is  baseu  on  a  model  by  Birkemeier 
and  Dalrymple  (19  76) . 

The  version  of  the  SSSP  most  suitable  for  use  on  the  HP- 
9845B  consists  of  a  main  program  which  controls  nine  subrou- 
tines.  A  flow  chart  of  the  program  is  given  in  Figure  1.   The 
subroutine  format  of  the  SSSP  allows  for  easier  test  and 
evaluation  by  component,  modification  or  program  correction, 
and  reduced  memory  requirements  on  the  mini  computer.   The 
main  disadvantage  of  the  subroutine  format  is  some  reduction 
in  calculation  speed. 

B.  REFRACTION  OF  SEA  SWELL 

The  wave  direction  is  determined  using  conservation  of 
wavenumber  for  a  steady  state  process  (Phillips,  1977): 

9(k  cos  9)    3(k  sin  9)      n       m 
9y      '      3x  U       {    ' 

where  (x,y)  are  the  on-offshore  and  alongshore  directions, 
and  9  is  the  counter  clockwise  angle  from  the  X  axis  to  the 
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Figure  1.   Sea  Surf  and  Swell 
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direction  of  swell  propagation  and  k  is  the  wavenumber .   It 
is  assumed  that  (1)  is  valid  in  both  open  water  and  within 
the  surf  zone. 

Noda  (19  72)  developed  an  iterative  scheme  to  solve  (1) 
for  the  wave  angle.   The  convention  for  grid  location  of  the 
numerical  model  is  shown  in  Figure  2 .   The  differential  equa- 
tion is  central-differenced  in  the  X-direction  and  forwarded- 
differenced  in  the  Y-direction.   The  finite  difference  form  of 
the  equation  usee  by  Noda  (19  72)  with  weighting  on  the  forward 
difference  term  added  to  increase  stability  (.after  Abbot, 
1979)  is  used  here: 


I.  .   =   sin"1^ — [(i)(k  sin  9).,,  .  . 

±  >  J 


-  (1  -  2t) (k  sin  0)   -  . 

(2 
+  (t) (k  sin  6) 


i+1, j+1 


^((k  cos  B)ifj+1  -  (k  cos  8)^^)]} 


where  weighting  factor  of  x  =  0.25  is  used.   The  initial  guess 
for  the  wave  angle  is  calculated  by  applying  Snell's  law  for 
parallel  bottom  contours.   (This  gives  an  exact  solution  in 
the  case  of  a  plane  sloping  beach.)   The  initial  guess  field 
of  9 .    is  iterated  until  9 .    is  within  a  half  percent  of 
the  previous  iteration. 
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Figure  2.   Differencing  scheme  used  in  the  model.   The 

grid  points  are  identified  by  subscripts  with 
the  first  subscript  indicating  position  from 
the  beach  with  increasing  numbers  farther 
offshore.   The  second  subscript  indicates 
position  alongshore  with  increasing  numbers  to 
the  right  when  the  beach  is  viewed  from  a 
ship  off  the  coast. 
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The  wave  numb  :  r  k.    is  the  positive  root  of  the  linear 
wave  dispersion  relationship 


a   =   (g  ktanh  kd)  1/2  (3) 


where  d  is  the  water  depth  and  a  is  the  angular  frequency. 
In  the  original  SSSP  model  by  Wang  and  Chen  (1983) ,  the  wave- 
number  was  solved  for  in  (3)  using  an  iterative  Newton's  method 
However,  to  speed  up  the  model,  the  wavenumber  is  calculated 
here  using  a  6th  degree  polynomial  fit  of  (3)  as  given  by 
Hunt  (1979) : 


k   =   -2—  (S  +  (1  +  0.66667S  +  0.35550S2 


/gd 


+  0.16084S3  +  0.06320S4  +  0.02174S5 


+  0.00654S6  +  0.00171S7  +  0.00039S8 


(4) 


+  0.00011S9)"1) "1 


where : 


S   =   a2(d/g) 


This  non-iterative  approximation  resulted  in  an  observed  90 
percent  savings  in  computation  time  on  the  HP-9845B  with  the 
same  accuracy. 
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C.   WAVE  HEIGHT  TRANSFORMATION 

The  information  of  wave  direction  is  used  to  calculate 
an  initial  guess  for  the  energy  field.   The  initial  energy  at 
each  grid  point  is  calculated  by  assuming  conservation  of 
energy  over  straight  and  parallel  contours. 


9(E  C   cos  9) 
2 

3x 


=   0  (5) 


Integrating  (5)  leads  to  the  initial  energy  field: 

C   cos  9 

g 

where  the  subscript  (o)  refers  to  conditions  at  the  offshore 
boundary  and  C   is  the  group  velocity.   The  initial  guess 
energy  field  is  used  to  calculate  the  wave  height  field. 

The  calculation  of  wave  height  is  based  on  the  energy 
balance  equation  including  refraction,  dissipation  and  growth. 
Three  assumptions  are  made  to  speed  calculation  and  reduce 
model  complexity  (Chen  and  Wang,  1983) .   First,  it  is  assumed 
that  wave  energy  within  a  particular  spectral  band  stays  in 
that  band,  which  allows  linear  superposition  of  wavelets.   The 
second  assumption  is  that  each  frequency  component  can  be 
described  by  a  single  mean  direction,  9~.   Third,  it  is  assumed 
that  there  are  no  wave-current  interaction.   With  these 
assumptions,  the  steady  state  energy  flux  equation  for  each 
spectral  band  can  be  written  (Longuet-Higgins ,  and  Stewart, 
1960,  1962) . 
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p  g  |tr(S(f)  C  cos  6)  +  p  g  f-.S(f)C  sin  0) 
ox.  g  ay       g 


=   -  ed(f)  +  £g(f)  (7) 


where: 

S(f)  =  wave  energy  density  at  spectral  frequency  f 

p  =  density  of  the  fluid 

g  =  acceleration  due  to  gravity 

e,  (f)  =  energy  dissipation  through  bottom  friction 

e  (f)  =  energy  generation  due  to  wind  stress 

The  energy  generation  term,  e  ,  is  based  on  the  mechanisms 
used  for  wind  energy  transferred  to  the  wave  field  described 
by  Phillips  (1957)  and  Miles  (1957 ; 1959a, b; 1962) .   The  Phillips 
mechanism  is  due  to  random  atmospheric  pressure  fluctuations 
advected  by  wind  resonantly  interacting  with  a  perturbed  sea 
surface.   The  wave  energy  generated  by  the  Phillips  wind 
mechanism  is 


S    =   p  ga  (8) 


where 


4tt  k  a, 
a   =   j-gi  P(k,a)  (9) 

p  g 


and 


a,   =   a/g  (10) 
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Barnett  (1968)  fitted  the  following  empirical  relationship 
for  P  based  on  field  measurements : 


P(k,o>   =   6.13  xio  w    2 

it  a      v2  +  (k  sin  6) 


vl 
[~2 ~ j]       (11) 

v,  +  (k  cos  6  -  A) 


where: 

6   =   angle  between  wind  and  wave 
w  =  wind  speed  (meters/second) 

A   =   a/w 

n  -x       1-28 

v,   =   0.3 

A  cro   0.95 

v~   =   0.52 

This  mechanism  is  important  for  wave  initiation  from  calm 
water  and  results  in  a  linear  increase  in  energy  with  time 
(Hasselman,  1960) . 

The  Miles  energy  transfer  mechanism  is  important  only 
after  waves  have  formed  (LaBlond  and  Mysak,  1978) .   When  wind 
blows  over  an  irregular  sea  surface,  flow  separation  can 
occur  in  the  lee  of  the  wave  crest  creating  pressure  differ- 
ences.  Momentum  is  transferred  to  the  waves  as  a  result  of 
these  pressure  forces.   The  instability  mechanism  of  Miles* 
results  in  exponential  wave  growth  (Barnett  (19  6  8) ;  Hasselmann 
(1960) )  . 
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S        =      pgB[S(f)]  (12) 

m 


where: 


Asa  rW  cos    6        R1  , . 

~2¥[ 5 8]  (13) 


s   =   ratio  of  air  to  water  density 
A  =   5 
B   =   .9 


The  wind  wave  generation  is  then  the  sum  of  (S  )  and  (S  ) : 

p        m 


e    =   S   +  S  (14) 

s      p    ra 


The  wave  growth  is  limited  to  a  fully  risen  sea  at  each 
spectral  band.   The  limits  on  development  are  based  on  Phillips 
(1957)  equilibrium  spectrum  modified  by  Kitaigorodskii  et  a_l 
(1975)  and  Thornton  (1977)  for  the  shallow  water  condition: 


S(f,5)   =   B  g2f"5H(6)r(Wn)  (15) 


where : 


H(6)   =   (JL)  (cos4(5))  (16) 

J7T 


accounts  for  angular  spreading, 


9  2W^Z(W    )  -1  n7v 

HWn)       =      E"2(W){1   +  2 —  " }  un 

n  sinh[2WZZ(W    ) 

n        n 
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and 


d  1/2 
Wn  =  a(-)  (18) 


where 


E(Wn) 


satisfies  the  condition 


Z  Tanh(W  Z)   =   1  (19) 


r (W  )  =  1  in  deep  water,  and  for  Wn  less  than  1, 


and 


r(wn)  =  \*2n  (20) 


S(ff5)   =   jqgdf'3H(6)  (21) 


_2 
where  q  is  a  coefficient  on  the  order  of  10   ;  a  value  of 

.0  73  is  used  in  Wang's  model. 

The  bottom  dissipation  term,  £,,  is  due  to  work  done  on 

the  bottom  by  the  bed  shear  stress.   The  bed  shear  stress  is 

proportional  to  velocity  squared  and  the  rate  of  energy 

dissipation  is  proportional  to  the  cube  of  the  velocity. 

Hasselman  and  Collins  (1968)  proposed  a  quasi-linearized 
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expression  for  a  frequency-de   ident  dissipation  function: 


ed(f)   =   P  cf|Ub|  ub(f)  (22) 


where  c_  is  a  bed  shear  stress  coefficient,  uu(f)  is  tne  wave 
induced  bottom  velocity  at  a  particular  spectral  component, 
and  U,  is  the  bottom  velocity  for  the  total  flow  field. 
Applying  linear  wave  theory  transfer  functions  relation  the 
energy  spectrum  of  the  surface  elevation  to  the  velocity 
field,  (22)  can  be  expressed  as: 


2  2 

e,(f)   =  pc.  [-5 — 2* ]S(f)U,  (23) 

a  cosh  (kd) 


where  the  term  in  the  brackets  is  the  spectral  transfer 
function.   The  total  bottom  velocity  is  obtained  by  applying 
the  spectral  transfer  function  to  the  entire  spectrum: 


oo  2    2  1/2 

(U.  )       =       [    /         J  k    %U)    df]  (24) 

D  0        a     cosh    kd 


D.   THE  SURF  MODEL 

This  model  describes  the  surf  zone  width,  number  of  lines 
of  breakers,  breaker  type,  breaker  angle,  significant  and 
maximum  breaker  height,  and  the  effective  surf.   The  calcula- 
tions are  based  on  empirical  formula  given  in  the  Shore  Pro- 
tection Manual  (U.S.  Corps  of  Engineers,  1977).   The  input 
wave  is  monochromatic,  and  requires  only  significant  wave 
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height  and  period.   Initial  wave  height  and  direction  are 
determined  for  the  entire  model  grid  using  swell  direction 
and  wave  height  calculations  described  in  Chapter  II,  sections 
B  and  C.   The  depth  of  the  first  breaker  is  determined  by 
comparing  the  calculated  wave  height  with  the  incipient  break- 
ing criteria  (Noda  et  al . ,  1974;  Thornton  and  Guza,  1983). 


Hb  dh 

-~  =   0.43  Tanh(-S-)  (25) 

J-i  L 


The  width  of  the  surf  zone,  X,  ,  is  then  calculated  using 
the  relationship  between  breaker  height,  breaker  depth,  wave 
period,  and  bottom  slope  (Weggel,  1972). 


Hb(x)  Hb(x) 

-,  ,  »   =   b(m)  -  a(m)  5—  (26) 

db(x)  ^T- 


assuming 


H,  (x)   =   c  +  dX,  (26a] 

b  b 


db(x)   =   e  +  fXb  (26b) 


Solve: 


,2 


(a'bf)Xj:  +  (a'ed  +  a'cf  -  b'f  +  d)  2^  +  (a'ce  +  c-b'e)   =   0 

(26c) 
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The  coefficients  c,  d,  e  and  f  are  determined  by  using 
standard  slope  intercept  techniques.   The  positive  root  of 
this  quadratic  equation  is  the  width  of  the  surf  zone. 


a(m)   =   1.36(1  -  e"19m)  (27) 


b(m)   -   I   ''-19.5m  (28> 

1  +  e 


a'   = 


a(m) 


b'   =   b(m) 

where : 

m  =  bottom  slope  between  grid  points 

The  first  breaking  wave  height,  denoted  as  HBK  in  the 
program,  is  then  determined  using  the  surf  zone  width  by 
interpolating  between  grids. 


HBK   =   ((HI  j((2-I)  +  (HI_1  jKl-D)  +  (  I,J  AxI"1/J)x 

(29) 


where  H  is  the  original  wave  height  from  energy  considerations 
only  and  I  and  1-1  are  the  grids  that  bracket  the  deepest 
breaker  in  the  on-offshore  direction. 

The  depth  of  the  first  b~=aker,  DBK,  is  calculated  also 
using  the  surf  zone  width  ar.   ihe  interpolation  between  grids: 
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D    -D 
DBK   =   [(D) (2-1)  +  (DI_1^J)  (1-1)]  +  (  lfJ    Ax   /J)X 

(30) 

where  D  is  the  bottom  depth  in  the  program. 

Additional  lines  of  breakers  shoreward  of  the  initial 
breaker  are  calculated  by: 


b'  d, 
H        =       b  next  mi 

"b  next     1  +  a'  d.     7  *JJJ 

b  next 


where 


db  next  =  hI,J  +  hl+1^x"hl-J[X-(I-l)Ax]      (32) 


The  angle  of  the  breaking  wave  (Angz)  is  calculated  by 
interpolating  the  distance  between  grid  points  and  the  linear 
change  in  the  wave  direction: 


Ang    -Ang 
Angz   =  An9l  j  +  ( ±^-^ ±*±)   (X  -  (1-1) Ax)       (33) 


where  Ang  is  the  calculated  wave  direction  at  the  grid  point. 
The  type  of  breaker  is  based  on  the  surf  parameter: 


£h   "   m  T75  (34 

b      (H  /L  )T72 
o  o 


where : 
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m  =   bet torn  slope 

H   =  deep  water  wave  height 

2 
L   =   1.56T  ,  deep  water  wave  length 


Breaker  type  is  specified  by  the  following  criteria 


e.   >   2  surging  or  collapsing 


0.4  <  e.   <   2  plunging 


e.   <   0.4  spilling 


Effective  surf  is  calculated  in  accordance  with 
COMNAVSURFPAC/COMNAVSURFLANT  Instruction  3840.1  for  use  in 
naval  amphibious  operations.   "Effective  surf  height  is  a 
numerical  value ,  expressed  in  feet,  and  provides  a  guide  for 
judging  the  feasibility  of  landing  operations  for  any  given 
landing  craft  under  existing  or  forecast  conditions."   Effec- 
tive Surf  is  a  function  of: 

1.  Breaker  Height 

2.  Breaker  Period 

3.  Breaker  Type 

4.  Littoral  Current 

5.  Wind  direction  and  velocity 

6.  Secondary  wave  height 

7.  Breaker  angle 

The  parameters  for  the  empirical  calculation  of  effective 
surf  are  calculated  in  the  surf  condition  part  of  the  model 
except  for  littoral  current.   The  nearshore  current  circulation 
in  the  model  is  described  in  the  next  section. 
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E.   NEARSHORE  CIRCULATION  MODEL 

The  nearshore  circulation  is  generated  by  changes  in  momen- 
tum due  to  breaking  waves  within  the  surf  zone  and  by  wind 
induced  surface  stress.   The  nearshore  circulation  model  is 
based  on  the  depth  integrated  and  time  averaged  equations  of 
motion  by  Phillips  (1977)  .   The  breaking  wave-induced  currents 
are  calculated  by  alternately  solving  the  momentum  and  the 
continuity  equations.   The  equations  are  solved  numerically 
with  the  finite  differences  scheme  taken  directly  from  Kirby 
and  Dalrymple  (1982)  which  is  an  updated  version  of  Birkermeier 
and  Dalrymple  (19  77) .   The  x  and  y  momentum  equations  are 
integrated  over  the  depth  and  time  averaged  over  one  or  several 
wave  periods.   The  equations  contain  bottom  shear  stress,  mean 
surface  shear  stress  and  excess  mean  momentum  stress  due  to 
wave  action.   To  simplify  the  model,  the  formulation  neglects 
the  nonlinear  mean  momentum  flux  terms  and  the  lateral  (hori- 
zontal) transfer  of  turbulent  momentum.   The  continuity  equation 
is  given  by: 


|£  +  |^[u(d  +  n)J  +  |y[vtd  +  n)J   =0  (35) 


where  r\    is  the  free  surface  elevation,  d  is  the  still  water 
depth,  (u,v)  are  the  (x,y)  horizontal  velocity  components  and 
the  overbar  indicates  time  averaging.   The  horizontal  momentum 
equations  are 

3S     3S 

3t   "   "  9  at  "  — — =r(-^x-+-9/-Tsx+Tbx)         (36) 

p(d+n)        2 
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go        gc 

»t        9  3y    p(d+_  (  3x      3y     xsy    Tby)      (37) 


The  S   ,  S   ,  S   ,  S    are  the  mean  momentum  flux  due  to  the 
xx   xy   yx    yy 

wave  motion  (commonly  referred  to  as  radiation  stresses) . 

t   ,  t    are  the  surface  shear  stresses  due  to  wind,  and 
sx   sy 

t,  ,  t,   are  the  bottom  stresses, 
bx'   by 

The  radiation  stresses  are  given  by  (Longuet-Higgins  and 
Stewart,  19  64) : 


S     =   E[(2n  -  l/2)cos29  +  (n  -  l/2)sin20]  (38) 


S    =   E[(2n  -  l/2)sin29  +  (n  -  l/2)cos29]  (39) 


S     =   S     =   |nsin(29)  (40) 

xy      yx     2 


1    2 

where  E  is  wave  energy  =  o-pgH  ;  n  =  C  /C. 

The  surface  wind  shear  stresses  are  calculated  (Van  Dorn, 
1953)  by: 


T     =   pK  W  W  (41) 

sx      K  '   '  X 


t     =   pK|w|W  (42) 

sy      K  '  '  y 


where  W  ,  W  are  the  wind  speed  in  the  (x,y)  directions,  and 
x   y 

—  6 
K  is  the  wind  stress  coefficient  of  order 
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The  bottom  shear  stresses  are  linearized  by  assuming  weak 
mean  currents  (U,V)  (LaBlond  and  Mysak,  19  78) . 


Tbx   "   p  W   cf  um  U  (43» 


Tby  ■   p  7  cf  um  V  t44» 


wh 


ere  cf  is  the  Darcy-Weisbach  bottom  friction  factor  and 


u   is  the  maximum  wave  orbital  velocity  at  the  bed: 

m  2 


H  a  ,,,. 

um     2  sinh  kd  *"■ 


The  maximum  time  step  At  is  based  on  the  linear  stability 
criterion  (Kirby  and  Dalrymple,  1982)  : 


At   <   n/((Ax)2  +  (Ay)2/2gd)  (46) 

However,  in  practice,  Kirby  and  Dalrymple  (19  82)  recommend 
using  a  quarter  of  this  value.   The  actual  At  used  in  the 
program  is 

2      2  1/2 

At   =  \    (Ax   +  $L  1 (47) 

(gd)x/z 

The  finite  difference  form  of  (36)  and  (37)  are  given 
by: 
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s      -s 
un^    ■  u^    +  Att<n_  ,  T-nT  T)^+  (  X  I"1,J  XXl-J; 

newi,J      oldI,J         X    lfJ         lfJ   AX  Axpd 


where : 


_  ,Sxyi,j+i  "Sxyi,j-i  +Sxyi-iyJ-n   Sxyi-i,j-i 

4Aypd  4Aypd         ' 

T         +  T  T,        +  T. 

SX  SXT-1  T        bXT  T      bXT-l  T 

+  '    I,J  2pd      '  >  -  <    —  2pd        "     l48: 


dI   J-l   +   dI"l   J-l 

d  =   ( — 2 '' )   is  t^e  average  depth  in 

the  X-direction 


v       =  v  ,  ,    +  At[(n,  ,  t  -  nT  ,)  S- 

new_  _      old_  _        I, J-l     I,  J   Ay 

X  z  J  J.  ,  J 


where : 


,Sxyi+i,j-i~Sxyi-i,j-i,  fSxyi-n,j"Sxyi-i,j. 

K  4Axpd         '~[  4Axpd      ' 


1       Aypd      '  +  l         2pd       ; 


byi  t      byi  T-i 

-   < 9nH L^-L)]  (49) 

2pd 


dI  J-l  +  dl'J"2 
d  =   (  lfJ   *  x ) 
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The  finite  difference  form  of  (35)  solving  for  n  is 


J.  /  J  J-  *  <-> 

-  n      (dI,J-l  +  dI+l,J-l,         ,dI,J-2+di,j-l 
UI+l,Jl        2Ax        '  +   I, J1      2"Ay" } 


dI  J  +  dT   T-l 

VI,J+1^       2Ay     n  ^UJ 


The  new  water  depth  is  then  given  by 


old   _new 
dI,J   =   (dI,J  "  ni,JJ     +  nI,J 


The  nearshore  circulation  solutions  are  "spun  up"  to 
steady-state  condition.   The  model  is  started  from  state  of 
rest  with  zero  wave  field  present  (Kirby  and  Dalrymple,  1982) 
To  reduce  the  effect  of  seiching  or  instability,  the  wave 
height  H  at  the  offshore  grid  is  gradually  brought  up  ~o 
its  full  value  using: 


H   =  Hq  tanh(^)  (51) 

s 


where : 


t   =  model  time 


T   =   arbitrary  fixed  period 

H    =   full  height  of  boundary  wave 
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The  hyperbolic  tangent  function  is  used  to  control  the 
growt.  of  the  boundary  wave  during  an  iterative  process.   The 
quantity  [2t/T  ]  increases  linearly  while  the  tanh  function 
increases  smoothly  toward  unity. 


31 


III.   MODEL  APPLICATION  ON  THE  HEWLETT-PACKARD  98  4  5B 

A.  CONVERSION  INTO  BASIC 

The  U.S.  Navy  contracted  for  the  development  of  the  wave 
and  surf  model  which  was  written  in  FORTRAN.   In  its  original 
form,  the  model  was  written  for  a  large  non-portable  computer. 
To  make  the  model  useful  tactically  for  mobile  U.S.  Navy  forces, 
the  program  is  rewritten  specifically  for  deployed  micro 
computers  like  the  HP-9845B.   It  is  required  that  the  model 
be  written  in  a  user  friendly  format  with  common  naval  conven- 
tions and  terms.   Error  checking  and  input  review  had  to  be 
programmed  into  the  model  to  be  operationally  acceptable  to 
the  U.S.  Navy.   The  Wave  and  Surf  Model  had  to  run  on  the 
HP-9845B  to  give  reliable,  useful  and  easily  understood  output. 

B.  MODIFICATIONS  IN  THE  WAVE  AND  SURF  MODEL  FOR  THE  HP-9  845B 
The  original  contractor  delivered  a  FORTRAN  version  of 

the  wave  and  surf  model  written  in  a  format  of  "stand  alone 
subroutines."   This  version  was  suitable  for  conversion  to 
the  HP-9845B  because  its  size  was  within  the  micro  computer's 
memory  requirements.   The  subroutine  format  also  allowed  for 
easier  testing  of  the  model.   The  main  change  made  in  the 
program  involved  the  use  of  common  statements.   The  FORTRAN 
version  made  extensive  use  of  labeled  common  blocks.   The  HP 
basic  only  allows  one  common  statement.   This  minor  problem 
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was  solved  by  increasing  the  size  of  variable   n  the  call 
statement  to  most  subroutines.   Only  the  most  ; ften  used 
general  variables  were  used  in  the  single  common  statement. 

The  language  of  the  input  statements  is  changed  to  agree 
with  naval  conventions.  Wind  and  swell  are  described  by  the 
direction  they  arrive  from.  The  beach  is  described  from  the 
point  of  view  of  a  ship  offshore. 

The  interactive  subroutines  are  expanded  with  lines  of 
code  that: 

1)  explain  the  input  needed 

2)  ask  for  the  data  input 

3)  error  check  the  data  input  before  accepting  or 
rejecting  and  requesting  correct  input. 

The  program  was  expanded  in  the  calculation  code  with 
explicitly  assigned  variables  for  counters  and  matrices.   For 
example,  effective  surf  was  stored  into  a  unique  array  rather 
than  sharing  an  array  used  for  several  other  calculations. 
Empirical  values  from  the  COMNAVSURFPAC/COMNAVSURFLANT  Instruc- 
tion 3840.1  were  put  into  an  array  rather  than  use  the  FORTRAN 
version  data  fitting  calculation.   The  additional  memory 
space  for  these  arrays  was  traded  for  explicitly  assigned 
data  and  clearer  programming. 

Labels  were  placed  before  every  significant  calculation 
in  the  model  code.   Spaces  were  added  between  code  to  separate 
functional  groups.   Indentations  of  the  code  were  made  to  show 
loops  and  double  loops.   These  changes  added  to  the  readability 
of  the  program  and  made  the  model  logic  easier  to  follow. 
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In  addition,  changes/corrections  were  made  to  the  sea 
surf  and  swell  model  for  undefined  constants,  incorrectly 
dimensioned  arrays,  logic  errors  and  improperly  defined 
angles.   Some  of  these  changes  were  made  with  the  cooperation 
of  the  contractor.   In  the  program  code,  each  change  is 
documented. 
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IV.       DISCUSSION 

A.  USABILITY 

The  sea,  surf  and  swell  model  is  programed  to  be  inter- 
active in  language  a  sailor  can  understand.   The  units  for 
wave  height  and  bottom  depth  entires  can  be  either  meters  or 
feet.   Grid  spacing  can  be  in  nautical  miles  only  for  the 
open  water  version.   However,  the  nearshore  model  allows 
either  feet  or  meters.   All  entries  are  requested  with  a 
format  example.   Where  possible,  the  format  example  is  in 
standard  World  Meteorological  Organization  code,  which  should 
be  familiar  to  all  U.S.  Navy  Aerographers  Mates.   Suggested 
input  values,  or  default  values,  are  offered  where  appropriate 
to  make  the  model  more  "user  friendly." 

B.  DATA  INPUT  REQUIREMENTS 

The  quality  of  the  input  data  is  the  most  important  factor 
determining  accurate  wave  and  surf  forecasts.   The  bathymetry 
begins  to  affect  sea  and  swell  when  the  depth  is  less  than  half 
a  wave  length.   For  the  same  depth  of  water,  the  effect  of 
the  bottom  will  be  strongest  on  longer  period  waves.   The 
more  shallow  the  water,  the  greater  the  effect  will  be  on  all 
waves.   Accurate  bathymetric  data  is  needed  in  deeper  sections 
of  the  model  grid,  but  becomes  increasingly  important  for 
shallow  waters.   The  wave  period  and  depth  is  used  to  compute 
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wavelength  and  phase  speed.   These  parameters  determine 
forecast  wave  direction,  energy,  surf  type  and  currents 
generated  in  the  surf  zone.   For  the  same  reason,  initial 
wave  direction  is  as  important  as  the  wind  direction.   The 
bottom  type  input  influences  energy  dissipation  and  surf  zone 
currents.   The  surf  model  calculation  is  least  sensitive  to 
the  bottom  type. 

The  most  difficult,  and  clearly  the  most  important  input, 
is  the  offshore  wave  input.   Careful  consideration  must  be 
given  to  the  many  sources  of  deep  water  waves.   The  open 
model  sums  the  energy  from  all  deep  water  wave  inputs.   The 
wave  spectrum  input  must  be  accurately  described.   The  ocean 
can  be  nearly  flat  concealing  long  period  swell  that  can 
generate  catastrophic  surf  on  beaches.   The  ocean  can  have 
fully  developed  seas  generated  by  continuous  wind.   In  a 
fully  developed  sea,  the  energy  of  the  waves  is  distributed 
over  a  wide  range  of  periods  that  must  be  entered  into  the 
model  for  accurate  surf  forecasts.   Only  with  a  sound  knowledge 
of  deep  water  wave  conditions  can  the  model  be  used  to  fore- 
cast shallow  water  conditions. 

The  open  water  part  of  the  model  has  the  most  options  for 
wave  energy  input,  which  include: 

1.  A  single  narrow  band  swell. 

2.  A  significant  wave  height  and  period  which  can  be 
expanded  in  terms  of  a  Bretschneider  spectrum. 

3.  Specific  energy  levels  can  be  entered  at  appropriate 
frequencies  such  as  from  the  Navy  Spectral  Ocean  Wave 
Model. 
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It  is  important  to  note  that  the  output  from  this  model  is 
a  single  significant  wave  heirnt  and  period  moving  in  the 
direction  of  maximum  energy.   This  could  be  a  poor  represen- 
tation if  the  seas  were  large  and  from  diverse  directions. 
The  surf  and  nearshore  circulation  parts  of  the  model 
allow  only  a  monochromatic  wave  input.   This  input  could  be 
directly  from  the  open  water  calculation  or  an  input  of  the 
user's  choice.   If  there  is  more  than  one  dominant  wave  group 
approaching  the  beach,  the  model  representation  of  the  surf 
zone  will  be  incomplete. 

C.   SELECTING  THE  GRID  SIZE 

The  grid  size  and  the  distance  between  grid  points  will 
determine  the  model's  ability  to  resolve  the  wave  climate. 
The  open  water  model  has  grid  distances  in  nautical  miles. 
The  grid  mesh  should  cover  the  area  from  deep  water  waves 
into  shallow  water,  outside  the  surf  zone.   The  grid  size  for 
the  surf  and  nearshore  is  measured  in  feet,  and  the  selection 
is  critical  to  resolving  the  surf  zone.   Kirby  and  Dalrymple 
(1982)  outline  the  following  requirements  for  the  nearshore 
circulation  grid. 

1.  The  grid  must  extend  offshore  far  enough  to  remove  the 
offshore  region  of  the  domain  from  the  influence  of 
currents  driven  by  the  surf  zone,  and  to  allow  for 
the  specification  of  a  uniform  longshore  depth  which 
will  not  significantly  alter  the  wave  refraction 
results  in  the  nearshore. 

2.  The  grid  mesh  must  be  fine  enough  to  resolve  the  surf 
zone  adequately. 

3.  In  the  event  that  the  effect  of  a  single  physical 
feature  isolated  in  the  longshore  direction  is  to  be 
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modeled,  the  longshore  extent  of  the  grid  must  be 
sufficiently  large  to  isolate  the  physical  system 
from  the  effect  of  images  created  by  the  longshore 
periodicity  requirement. 

Experience  running  the  model  suggests  that  a  minimum  of 
3  grid  points  must  be  inside  the  surf  zone.  This  would  not 
be  possible  for  waves  breaking  directly  on  the  beach.  A 
further  constraint  is  the  number  of  grids  in  this  model  (10 
offshore  rows  and  13  columns  alongshore)  due  to  the  limita- 
tions of  speed  and  storage  area  of  the  micro  computer. 
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V.   MODEL  TEST  AND  EVALUATION 

A.   DEEP  WATER  WIND  GENERATED  WAVES  COMPARED  WITH  JONSWAP 
CURVES 

The  model  approximates  the  area  under  the  spectral  density 
curve  (wave  energy)  with  rectangles.   The  more  rectangles 
used,  the  better  the  approximation  will  be,  but  computational 
time  also  increases.   The  open  water  model  requires  approxi- 
mately 4  5  minutes  to  calculate  each  frequency  band  on  the 
HP-9845B.   Narrow  band  swell  waves  might  be  closely  approximated 
with  a  single  energy  band  and  a  band  width  of  .01  Hz.   A 
broad  band  wave  spectrum,  such  as  wind  waves,  needs  to  be 
represented  by  multiple  energy  bands. 

To  test  the  wind  generation  model  and  the  effects  of 
bottom  friction  and  band  width  selection,  the  model  was  run 
for  constant  depths  of  100,  10,  and  5  meters  (Figs.  3,  4,  and 
5).   The  initial  wave  height  is  zero.   The  bottom  friction 
coefficient  chosen  is  .01  for  a  fine  sand  bottom.   The  grid 
spacing  is  10  nautical  miles  resulting  in  a  maximum  fetch  of 
90  nautical  miles.   The  results   are  compared  with  the  JONSWAP 
significant  wave  height  as  a  function  of  fetch  from  Bishop 
(1983).   The  model  results  favorably  agree  with  the  JONSWAP 
curve  considering  the  relative  coarseness  of  the  energy 
bands.   For  the  100  m  depths,  bed  friction  was  insignificant 
in  changing  the  wave  height  as  expected. 
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Figure  3.   Significant  wave  height  generated  by  a  15  kt 
wind  in  100  m  depth.   Model  computed  values 
(solid  line)  are  the  sum  of  spectral  contri- 
butions at  .05,  .10,  .15  and  .20  Hz  frequency 
bands.   The  dashed  curve  represents  wave 
height  generated  by  the  JONSWAP  spectrum  as  a 
function  of  fetch  and  wind  curves  in  deep  water 
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Fig.  4.   Significant  wave  height  generated  by  a  15  kt 

wing  in  10  meter  depth.   Model  computed  values 
(solid  line)  are  the  sum  of  spectral  contribu- 
tions at  .05,  .10,  .15  and  .20  Hz  frequency 
bands.   The  dashed  curve  represents  wave  height 
generated  by  JONSWAP  spectra  as  a  function  of 
fetch  and  wind  in  deep  water. 
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Figure  5.   Significant  wave  height  generated  by  a  15  kt 
wind  in  5  meter  depth.   Model  computed  values 
(solid  line)  are  the  sum  of  spectral  contribu- 
tions at  .05,  .10,  .15  and  .20  Hz  frequency 
bands.   The  dashed  curve  represents  wave 
height  generated  by  JONSWAP  spectra  as  a 
function  of  fetch  and  wind  in  deep  water. 
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The  model  was  run  for  the  same  input  with  the  flat  bottom 
depth  changed  to  10  meters.   The  effect  of  friction  only 
slightly  reduced  the  fully  arisen  significant  wave  heights. 
However,  the  model  generated  fetch  limited  wave  heights  were 
approximately  half  a  foot  higher  than  the  JONSWAP  results. 
The  increase  in  the  height  of  the  "shallow  water  waves"  was 
probably  a  result  of  the  change  to  a  shallow  water  saturation 
spectrum.   In  shallow  water  the  slope  of  the  energy  saturation 
spectrum  with  respect  to  frequency  is  -3  (Thornton,  19  77) . 
In  deep  water,  the  slope  is  -5  (Phillips,  1966).   This  results 
in  relatively  more  wave  energy  at  higher  frequencies  for  the 
shallow  water  spectra. 

To  test  the  bottom  friction  in  the  model,  a  run  was  made 
with  the  same  input  as  before  except  the  depth  was  changed 
to  5  meters.   The  significant  wave  height  curve  shows  a  slow 
decrease  after  the  first  grid  point.   The  model  calculates 
larger  frictional  dissipation  as  the  down  grid  distance  in- 
creases.  The  model  nearly  reached  a  constant  wave  height, 
balancing  generation  and  dissipation  forces  after  80  nautical 
miles.   The  values  of  predicted  wave  heights  by  the  JONSWAP 
spectrum  agree  with  the  present  model  within  10  nautical  miles, 
but  tend  to  be  larger  for  a  longer  fetch.   For  an  open  sea  in 
a  shallow  region,  the  JONSWAP  spectrum  ignoring  bottom  friction 
effect  gives  a  wave  height  four  times  larger  than  the  pre- 
dicted by  the  present  model;  this  result  emphasizes  the 
importance  of  including  bottom  friction  in  shallow  water  wave 
models . 
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The  wind  generation  part  of  the  model  was  again  tested 
for  30  kt  winds  blowing  over  the  same  area.   The  significant 
wave  heights  computed  by  the  model  increase  rapidly  as  the 
downwind  fetch  is  longer.   The  model  predicts  a  near  fully 
arisen  sea  after  only  50  nautical  miles.   The  JONSWAP  model 
predicts  approximately  linearly  increasing  wave  height  values 
for  30  kt  winds  until  reaching  fully  arisen  seas  after  200 
nautical  miles  of  fetch.   The  large  difference  between  the 
model  calculations  at  30  kts  and  the  JONSWAP  curves  suggest 
the  model  over-builds  seas  for  this  wind  speed  and  grid  spacing 

The  energy  transport  calculation  in  the  model  was  found 
dependent  on  grid  spacing  and  is  invalid  for  large  grid 
spaces  and  higher  frequencies.   Negative  energies  resulted  for 
input  at  frequencies  greater  than  .3  Hz  with  a  grid  spacing 
of  10  nautical  miles.   By  decreasing  Ax  =  Ay  to  one  nautical 
mile,  the  model  can  compute  energy  transport  for  .4  Hz  and 
higher.   A  "fix  up"  in  the  program  is  made  by  setting  any 
negative  energy  values  to  zero.   The  inability  of  the  model 
to  compute  energy  transformation  for  the  short  waves  should 
be  noted,  even  though  this  is  the  low  energy  side  of  the 
spectra. 

B.   NEARSHORE  CIRCULATION  ON  A  SLOPING  BEACH 

The  nearshore  current  prediction  part  of  the  model  was 
tested  on  a  uniform  sloping  beach.   Since  the  depth  is  uni- 
form in  the  alongshore  (Y)  direction,  the  wave  height,  wave 
direction,  current,  wave  set-up  and  surf  should  also  be 
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Figure  6.   Significant  wave  height  generated  by  a  30  kt 

wind  in  100  meter  water  depth.   Model  computed 
.values  (solid  line)  are  the  sum  of  spectral 
contributions  at  .05,  .06,  .07,  .08,  .09,  .10, 
.11,  .12,  .13  and  .14  Hz.   The  dashed  curve 
represents  wave  height  generated  by  JONSWAP 
spectra  as  a  function  of  fetch  and  wind  in 
deep  water. 
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uniform  in  the  Y  direction.   The  calculated  wave  directions 
were  constant  in  the  Y  direction,  the  wave  height  and  surf 
changed  a  few  percent  but  the  variations  of  current  and  wave 
set-up  in  the  Y  direction  deviate  more  than  20  percent. 
The  wave  height  and  surf  were  slightly  inconsistent  with  ex- 
pected results  and  this  is  thought  to  be  caused  by  the  finite 
differencing  scheme  producing  errors  from  the  boundaries.   The 
larger  variation  in  the  current  and  wave  set-up  values  is 
probably  due  to  the  inconsistent  mass  flux  calculations.   The 
boundary  errors  induced  in  the  wave  calculations  resulted  in 
a  nonuniform  forcing  term,  and,  therefore,  they  may  not  be 
able  to  describe  a  constant  current.   The  model  is  expected  to 
give  best  results  in  the  interior  of  the  computational  domain 
for  wave  and  surf  prediction. 

The  nearshore  current  calculation  portion  of  the  model 
was  compared  with  data  acquired  at  Santa  Barbara  on  4  February 
1980.   Initially  the  currents  calculated  were  one  order  of 
magnitude  smaller  than  the  observations. 

Consequently,  three  basic  problems  with  the  surf  model 
were  identified: 

1.  The  model  computed  only  36  seconds  of  surf  after 
beginning  with  still  water.   The  contractor's  program 
calculated  only  10  time  steps,  or  3.6  seconds.   This 
length  of  time  was  too  short  to  even  approach  a 
steady  state  current.   Kirby  and  Dalrymple  (1982) 
always  ran  their  nearly  identical  version  of  this 
model  for  a  1000  time  steps  representing  more  than 

5  minutes  real  time. 

2.  The  resolution  was  too  coarse  to  describe  the  current 
field.   The  finite  difference  scheme  is  unable  to 
compute  currents  nearer  than  2  gird  spaces  from  a 


46 


depth  of  zero  (at  the  beach  presumably) .   The  surf 
zone  for  this  test  was  closer  to  the  beach  than  2Ax 
or  20  meters . 

3.   The  changes  in  radiation  stress,  which  forces  the 

current,  were  very  small  across  the  breaker  line  where 
it  should  be  the  largest.   The  current  would  remain 
too  small  even  with  items  1  and  2  corrected  due  to 
weak  forcing.   Realistic  radiation  stress  was  not 
being  calculated  because  the  model  gave  incorrect 
wave  heights  for  breaking  in  the  surf  zone.   This 
correction  was  made. 

The  model  was  run  again  for  comparison  with  Santa  Barbara 

data.   The  grid  size  was  reduced  to  Ax  =  2.5  meters  and 

Ay  =  5  m.  The  time  step  for  this  small  grid  was  approximately 

.6  seconds.   The  HP-9845B  was  allowed  to  calculate  for  381 

time  steps  which  took  9  hours.   This  calculation  represented 

5.2  8  minutes  real  time  for  waves  entering  still  water.   The 

maximum  stable  time  step  is  a  function  of  the  grid  size 


TMAX   =   1/3 


gd 


The  current  was  approaching  a  steady  state  with  a  very  slow 
rate  of  convergence  (Fig.  7) .   It  would  take  more  than  23 
hours  of  CPU  time  on  the  HP-9845B  to  run  the  model  for  a 
1000  steps  as  done  by  Kirby  and  Dalrymple  (19  83) .   Attempts 
to  increase  the  time  step  and  decrease  the  CPU  time  resulted 
in  numerical  instability  (Figure  8) .   The  model  was  run  using 
the  same  input  except  the  At  was  increased  5  percent  each 
time  step  until  reaching  a  maximum  of  one  second.   Instability 
began  to  occur  after  20  seconds  with  values  oscillating  around 
the  stable  solutions.   After  25  seconds  the  unstable  model 
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Figure  7.  Current  generated  by  the  surf  as  predicted 
by  the  model  as  a  function  of  "spin  up" 
time.  A  total  of  3  81  time  steps  were  made 
representing  316  seconds  of  forcing  by  the 
surf.  The  dashed  line  is  an  extrapolation 
of  the  model's  calculation. 
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Figure  8.   Surf  zone  current  modeled  for  4  Feb  1980 

off  Santa  Barbara,  California.   The  stable 
increase  in  current  is  shown  with  the 
solid  line.   The  unstable  calculation  is 
shown  with  the  heavy  dashed  line. 


49 


of  the  current  became  large  and  changed  sign.   The  surf  model 
is  very  slow  to  converge  to  a  solution  and  unstable  for  even 
slightly  longer  time  steps. 

C.   SURF  CONDITION 

The  surf  zone  wave  heights  and  directions  calculated  by 
the  model  were  compared  with  observations.   The  comparisons 
were  made  with  observations  for  two  California  beaches  at 
Santa  Barbara  and  Torrey  Pines.   These  beaches  are  nearly 
plain  sloping.   The  offshore  boundary  condition  of  wave 
height,  period  and  direction  are  given  in  Table  I.   The 
beach  material  for  the  model  friction  coefficient  selection 
was  fine  sand  (Cf  =  .01).   The  modeled  and  observed  wave 
direction  agreed,  within  round-off,  to  the  nearest  degree 
for  all  days.   The  waves  turned  toward  normal  while  crossing 
the  surf  zone  exactly  as  expected. 

Model  comparisons  for  observations  on  3  February  are  shown 
in  Figure  (9)  which  shows  the  model  slightly  under-predicts 
the  wave  heights  but  follows  the  general  change  across  the 
surf  zone.   The  breaker  line  was  forecast  to  be  14.9  meters 
from  the  shore  with  a  height  of  .68  meters.   The  observed 
breakers  began  more  than  40  meters  from  the  beach  with  a 
height  of  .68  meters.   The  model  incorrectly  forecast  one 
line  of  breakers  when  at  least  two  were  observed. 

The  4  February  results  (Figure  (10))  show  reasonable 
prediction  of  wave  height.   But  using  the  wave  height 
information,  the  model  did  not  predict  the  location  or  height 
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Wave   direction 

(m) 

(m) 

(sec) 

C°) 

TABLE  I 
Offshore  Boundary  Conditions 

Date 

Santa  Barbara,  Leadbetter  Beach,  California 

3  Feb  80   4.1  .55        14.3          -7.8 

4  Feb  80   4.2  .56       14.3         -9.0 

5  Feb  80   4.1  .45        12.8         -8.4 

Torrey  Pines  Beach,  California 

4  Nov  78   6.7       .35        14.3  0 

10  Nov  78   6.6       .68        15.9  0 
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Figure  9.   Santa  Barbara  surf  zone  wave  height  on  3 
February  1980.   Model  computed  values  are 
represented  by  the  solid  line,  observed 
heights  are  dashed. 
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Figure  10. 


Santa  Barbara  surf  zone  wave  height  on  4 
February  1980.   Model  computed  values  are 
represented  by  the  solid  line,  observed 
heights  are  dashed. 
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or  breaker.   The  observation  indicated  two  lines  of  breakers, 
the  first  near  60  meters  and  the  second  near  38  meters.   The 
maximum  measured  wave  height  was  .67  meters  occurring  60 
meters  from  the  beach.   The  model  forecast  only  one  line  of 
breakers  .84  meters  high  and  16  meters  from  the  beach. 

The  5  February  results  (Fig.  (11))  show  the  model  almost 
forecast  the  surf  height  exactly  at  grid  points.   The  model 
forecast  one  line  of  breakers  11  meters  from  the  shore  with 
a  height  of  .61  meters.   The  observed  breakers  began  about 
39  meters  from  shore  with  a  height  of  nearly  .60  meters. 

The  model  was  also  compared  with  surf  zone  observations 
from  Torrey  Pines  State  Beach,  California  for  4  and  10 
November   1978  (Table  I) .   The  average  angle  was  near  normal 
to  the  shore.   Fig.  (12)  shows  the  model  input  boundary  wave 
on  4  November  height  was  too  small  by  about  .05  meters  (2 
inches) .   This  small  error  was  due  to  interpolation  between 
sensor  locations.   The  model  followed  closely  the  observed 
wave  heights  between  160  and  240  meters,  but  it  over-built 
the  wave  height  inside  160  meters.   The  model  placed  the  line 
of  breakers  82  meters  from  shore  with  a  height  of  .76  meters. 
Observed  breakers  were  near  160  meters  with  a  height  of  .48 
meters . 

The  10  November  1978  model  of  Torrey  Pines  was  based  on 
boundary  conditions  entered  at  6.6  meters.   Figure  (13)  shows 
uniform  over-prediction  of  the  wave  height  at  grid  points. 
The  model  forecast  two  lines  of  breakers,  the  first  line  133 
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Figure  11. 


Santa  Barbara  surf  zone  wave  height  on  5 
February  1980.   Model  computed  values  are 
represented  by  the  solid  line,  observed 
heights  are  dashed. 
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Figure  12. 


Torrey  Pines  surf  zone  wave  height  on  4 
November  1979.   Model  computed  values  are 
represented  by  the  solid  line,  observed 
heights  are  dashed. 
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Figure  13. 


Torrey  Pines  surf  zone  wave  height  on  10 
November  1979.   Model  computed  values  are 
represented  by  the  solid  line,  observed 
heights  are  dashed. 
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meters  from  the  beach  with  a  height  of  1.7  meters  and  the 
second  83  meters  from  the  beach  with  a  height  of  .75  meters. 
However,  the  observations  showed  no  corresponding  lines  of 
breakers.   The  observed  heights  of  the  lines  of  breakers  are 
.72  meters  or  less. 

The  wave  height  calculation  by  the  model  is,  in  general, 
close  to  the  observed  wave  height.   The  surf  zone  width 
based  on  (26)  is  incorrect  in  every  test  run.   The  resulting 
calculations  of  breaker  position,  angle  and  height  are 
dependent  on  the  width  of  the  surf  zone  and  therefore  these 
forecasts  are  incorrect. 
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VI.   SUMMARY  AND  CONCLUSIONS 

A.   SUMMARY 

A  numerical  model  describing  surface  gravity  wave  trans- 
formation in  shallow  water  region  is  written  in  BASIC  for  use 
on  the  HP-9845B.   The  original  version  of  this  model  was 
based  on  the  work  by  Wang  and  Chen  (19  83) .   The  model  consists 
of  an  open  water  wave  program  and  a  surf  program.   The  open 
water  model  considers  a  wave  energy  spectrum  generated  by 
winds;  the  surf  model  predicts  a  shallow  water  wave  field  and 
the  shear  currents  induced  by  wave  breaking  on  an  arbitrary 
bottom  topography.   The  model  is  tested  for  cases  of  constant 
depth,  a  uniformly  sloping  beach  and  actual  beaches.   The 
application  of  the  model  to  a  field  beach  indicates  the  poten- 
tial value  for  naval  operational  planning.   Combinations  of 
different  conditions  of  wind,  wave,  and  bottom  bathy  atry  can 
be  effectively  simulated  and  serve  as  the  basis  for  determin- 
ing an  amphibious  landing. 

The  open  water  model  is  used  to  forecast  sea  heights  over 
continental  shelf  and  wide  stretches  of  shallow  water  coast. 
The  calculations  of  wave  energy  are  the  result  of  wind  forcing, 
fetch  distance,  and  bottom  friction.   The  model  predicts  the 
spectral  distribution  of  wave  energy. 

The  surf  model  can  predict  wave  shoaling,  coastal  cur- 
rents and  multiple  breaker  lines.   In  the  model,  the  wave 
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heights  and  directions  are  calculated  for  a  section  of  beach. 
The  breaker  heights,  breaker  types,  and  location  of  severe 
breakers  are  empirically  determined  based  on  the  wave  height 
calculations.   This  provides  a  map  of  the  surf  zone  in  two 
dimensions  for  Navy  operations. 

One  of  the  drawbacks  of  this  model  is  that  the  model  is 
not  programmed  in  an  efficient  manner.   The  wavenumber  at  a 
local  water  depth  is  calculated  over  a  hundred  times  in  nearly 
every  subroutine  for  different  frequency  components.   To  allow 
for  various  options  chosen  and  setting  up  the  connectivity 
between  subprograms,  the  integrated  model  was  too  large  and 
not  very  well  structured.   Moreover,  the  unsteady  approach 
to  the  steady  state  solution  is  time  consuming  and  needs 
special  care  to  avoid  numerical  instability.   Using  a  steady 
state  approach  and  reforming  the  program  structure  would 
facilitate  the  application  of  this  model  and  improve  the 
operational  efficiency.   Further  modeling  efforts  should  be 
directed  toward  this  direction. 

A  problem  in  running  the  model  is  the  choice  of  the 
energy  bandwidth.   The  open  water  model  uses  energy  spectra 
to  describe  the  wave  field.   The  selection  of  energy  bands 
and  bandwidths  requires  knowledge  of  wave  spectra  and  gravity 
wave  theory.   The  alternative  to  user  choice  is  to  have  the 
model  calculate  energy  contribution  over  the  entire  wave 
spectra.   This  would  ensure  all  energy  contributions  across 
the  spectra  are  considered.   However,  CPU  time  would  increase 
to  more  than  6  hours  if  only  10  energy  bands  were  used.   (The 
Navy  Spectral  Ocean  Wave  Model  uses  15  energy  bands.) 
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B.   CONCLUSIONS 

This  sea,  surf  and  swell  model  at  present  is  not  suitable 
for  operational  use.   The  separate  regimes,  open  water  and 
nearshore  current  and  surf,  all  have  problems. 

The  open  water  part  of  the  model  requires  sophisticated 
spectral  energy  input.   The  average  Aerographer  Mate  does  not 
have  the  knowledge  to  properly  select  input  parameters .   Only 
when  the  sea  is  dominated  by  a  single  narrow  band  swell  with 
light  wind  is  the  wave  energy  input  fairly  easy  (input  wave 
height,  direction  and  period) .   In  addition,  the  open  water 
model  was  shown  to  build  up  the  sea  height  too  fast  for  wind 
speeds  of  30  kts  raising  doubt  concerning  the  wind  energy 
transfer  part  of  the  model. 

The  nearshore  current  part  of  the  model  is  very  slow  to 
converge  requiring  more  than  23  hours  to  complete  the  neces- 
sary 1000  iterations.   The  numerical  scheme  is  therefore 
unsatisfactory  for  the  HP-9845B  computer.   Recommendations 
include:   1)  using  a  more  stable  scheme  to  increase  the  time 
step  and  reduce  the  number  of  iterations;  2)  incorporate  an 
implicit  numerical  scheme  to  calculate  a  direct  solution  of 
the  currents  such  as  in  Wu  and  Liu  (1985) ;  or  3)  use  a  simple 
one-dimensional  current  model  implying  a  uniform  beach 
alongshore. 

The  surf  part  of  the  model  is  probably  the  most  "user 
friendly"  and  has  the  most  potential.   The  surf  part  of  the 
model  is  fast  running  and  gives  good  solutions  for  grid  point 
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wave  heights.   However,  the  width  of  the  surf  zone  calculation 
is  very  poorly  predicted.   This  error  in  turn  causes  the 
breaker  line  position,  height  and  angle  to  be  incorrect. 
The  incorrect  breaker  predictions  cause  the  effective  surf 
calculation  also  to  be  invalid. 
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